**************************************************************
* PROGRAM TO RUN AN INDIVIDUAL SET OF OLS REGRESSION RESULTS *
**************************************************************

capture program drop directctrlreg
program define directctrlreg

	set more off

	syntax [if/] [aweight], outcome(varname) [ctrlvars(namelist)] trimlevel(numlist max=1 >=0 <100) [header fh(name) printtoscreen]
	
	* Create list of control variables	
	local ctrlfvarlist ""
	foreach var of local ctrlvars {
		forval lag = 1/`lags' {
			if "`var'" == "survives" {
				local ctrlfvarlist "`ctrlfvarlist' i.`var'_l`lag'"
			}
			else local ctrlfvarlist "`ctrlfvarlist' c.`var'_l`lag'"
		}
	}

	* Work out if condition
	if "`if'" == "" local if "if 1"
	else local if "if (`if')"
	local if "`if' & (audityear >= 1999+4) & (yrsfiling >= 2)"

	* Work out trimming
	if (`trimlevel' == 0 ) {
		gen byte trimmed = 0
	}
	else {	//if trimlevel>0, construct variable noting which obs are trimmed
		* Cut (exclude) above a given percentile
		* Calculate percentile to trim at
		qui gquantiles `outcome' if (`outcome' > 0 & `outcome' < .), _pctile percentiles(`=100 - `trimlevel'')
		
		* Create trimmed variable
		gegen byte trimmed = max(`outcome' > `r(r1)'), by(utr_no audityear)
		
	}
	
	* Modify if condition to avoid having the same variable on both sides of the equals sign
	if ((inlist("`ctrlvars'","all","own")) | (("`ctrlvars'" == "survives") & ("`outcome'" == "survives"))) local if "`if' & (!inrange(yrssince,-4,-1))"
	* Modify if condition for trimming
	local if "`if' & (trimmed == 0)"

	* Create differenced outcome variable
	gegen double `outcome'_avg4 = total(`outcome'*inrange(yrssince,-4,-1)/4), by(utr_no audityear)
	qui replace `outcome'_avg4 = . if minyrssince > -4
	qui gen double `outcome'_diff4 = `outcome' - `outcome'_avg4
	drop `outcome'_avg4
	local outcomename "`outcome'_diff4"
	
	* Run OLS
	di `"qui reg `outcomename' i.yrssince20 i.yrssince20#i.treatment i.tax_year i.yrssince20#(`ctrlfvarlist') `if'  [`weight' `exp'], vce(cluster utr_no)"'
	qui reg `outcomename' i.yrssince20 i.yrssince20#i.treatment i.tax_year i.yrssince20#(`ctrlfvarlist') `if' [`weight' `exp'], vce(cluster utr_no)
	
	* Store results
	matrix eb = e(b)
	matrix eV = e(V)
	
	* Run margins
	qui margins, over(treatment yrssince20)
	matrix marginsb = r(b)
	matrix marginsV = r(V)
	
	*make it print nicely
	if (("`ctrlvars'" == "") local ctrlvars "NA"
	
	* Print header if requested
	if ("`header'" == "header") {
		if ("`printtoscreen'" == "printtoscreen") di as text %12s "method" %24s "outcome" %12s "ctrlvars" %12s "trimlevel" %10s "yrssince" %14s "ctrlgrp" %14s "ctrlgrpse" %14s "treatgrp" %14s "treatgrpse"  %14s "treateff" %14s "treateffse" %14s "pvalue" %14s "N"
		if ("`fh'" != "") file write `fh' "method,outcome,ctrlvars,trimlevel,yrssince,ctrlgrp,ctrlgrpse,treatgrp,treatgrpse,treateff,treateffse,pvalue,N" _n
	}
	

	* Print results
	forvalues yr20 = 15/28 {
		local ebcolnumb = colnumb(eb,"`yr20'.yrssince20#1.treatment")
		local marginsbcolnumbctrl = colnumb(marginsb,"0bn.treatment#`yr20'.yrssince20")
		local marginsbcolnumbtreat = colnumb(marginsb,"1.treatment#`yr20'.yrssince20")
		if ("`printtoscreen'" == "printtoscreen") {
			# delimit ;
			di as text %12s "ols" %24s "`outcome'" %12s "`ctrlvars'" %12.0g `trimlevel' %10.0g `yr20'-20
				%14.3f marginsb[1, `marginsbcolnumbctrl'] %14.3f sqrt(marginsV[`marginsbcolnumbctrl', `marginsbcolnumbctrl'])
				%14.3f marginsb[1, `marginsbcolnumbtreat'] %14.3f sqrt(marginsV[`marginsbcolnumbtreat', `marginsbcolnumbtreat'])
				%14.3f eb[1, `ebcolnumb'] %14.3f sqrt(eV[`ebcolnumb',`ebcolnumb'])
				%14.3f 2*ttail(e(N_clust), abs(eb[1, `ebcolnumb']/sqrt(eV[`ebcolnumb',`ebcolnumb'])))
				%14.0g e(N)
				;
			# delimit cr
		}
		# delimit ;
		if ("`fh'" != "") file write `fh' "ols,`outcome',`ctrlvars'," (`trimlevel') "," (`yr20'-20) ","
			(marginsb[1, `marginsbcolnumbctrl']) "," (sqrt(marginsV[`marginsbcolnumbctrl', `marginsbcolnumbctrl'])) ","
			(marginsb[1, `marginsbcolnumbtreat']) "," (sqrt(marginsV[`marginsbcolnumbtreat', `marginsbcolnumbtreat'])) ","
			(eb[1, `ebcolnumb']) "," (sqrt(eV[`ebcolnumb',`ebcolnumb'])) ","
			(2*ttail(e(N_clust), abs(eb[1, `ebcolnumb']/sqrt(eV[`ebcolnumb',`ebcolnumb'])))) ","
			(e(N)) _n
			;
		# delimit cr
	}
	matrix drop eb eV marginsb marginsV
	
	drop trimmed

	drop `outcomename'

end
